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ABSTRACT 


This  report  summarizes  the  work  done  by  the  SDL  during  the 
period  January  through  March  1968,  and  is  primarily  concerned 
with  seismic  research  activities  related  to  the  detection  and 
identification  of  nuclear  explosions  and  earthquake  phenomenon. 
Also  discussed  are  the  support  tasks  and  data  services  performed 
for  other  participants  in  the  VELA-Uniform  project. 
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Seismic  Data  from  Three  Nuclear  Explosions 


I.  INTRODUCTION 


This  quarterly  report  summarizes  the  technical  work,  support 
effort,  and  service  tasks  completed  during  the  period  January 
through  March  1968.  Current  or  past  work  is  mentioned  only  if 
it  relates  to  the  present  discussions. 

Reviews  of  technical  reports  completed  during  the  reporting 
period  are  contained  in  Section  II  under  descriptive  headings. 
Section  III  is  a  summary  of  the  support  and  service  tasks  per¬ 
formed  for  in-house  projects  and  for  other  VELA-Uniform  partici¬ 
pants.  This  report  concludes  with  Appendix  A  in  which  are  listed 
the  organizations  who  received  SDL  data  services  during  the 
period. 

II.  WORK  COMPLETED 

A.  LAS A  Travel-time  Anomalies  for  65  Regions  Computed  with 
the  Herrin  Travel-tine  Table,  Nov.  1966  Version 

A  report  was  published  concerning  the  results  of  computing 
P— wave  relative  travel— time  anomalies  at  all  subarray  center  in¬ 
struments  at  the  Large  Aperture  Seismic  Array  (LASA)  in  Montana. 

The  travel- time  table  used  for  computing  expected  times  is  the 

Herrin,  November  1966  version. 

A  total  of  626  teleseisms  are  used.  As  the  anomalies  are 
dependent  upon  event  distance  and  azimuth,  all  of  the  results  are 
separated  into  groups  such  that  the  subarrays  have  anomalies  which 
are  consistent  within  each  defined  region. 

The  following  description  explains  the  presentation  of 
the  results  in  Figure  1  using  the  reference  numbers  below. 
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Figure  1.  RELATIVE  TRAVEL-TIME  ANOMALIES 


1.  Source  of  expected  travel  timus.  In  this  report,  the 
Herrin  table,  November  1966  version,  is  used; 

2.  Reference  subarray,  R,  selected  for  computing  rela¬ 
tive  anomalies.  In  this  report,  all  anomalies  are  relative  to 
subarray  ^0.  The  following  relation  may  be  used  to  change  refer¬ 
ence  stations; 

A,  ,,  —  A.,  —  A,, 

i/j  i/r  j/r 


where  A  ,  is  the  anomaly  at  station  i  relative  to  a  new  reference 
i/j 

station  j. 

3.  All  expected  travel-times  in  this  report  have  been 
corrected  for  the  ellipticity  of  the  earth  such  that  the  computed 
anomalies  may  be  used  in  conjunction  with  other  programs  requir¬ 
ing  these  corrections. 

4.  An  arbitrary  geographic  name  given  to  the  event 

region. 

5.  Range  of  epicentral  distance  in  the  event  region. 

6.  Range  of  epicentral  azimuth  in  the  event  region. 

7.  Date  and  arbitrary  name  given  to  each  event. 

8.  Epicentral  distance,  in  kilometers,  from  the 
reference  subarray,  R. 

9.  Epicentral  azimuth,  in  degrees  measured  from  north 
to  east,  from  the  reference  subarray,  R. 

10.  Subarray  designator,  i. 

11.  Measured  travel-time  anomaly,  in  seconds,  at  sub¬ 


array  i  relative  to  subarray  R  for  the  kth  event; 

Ak,  =  Tk  -  Tk  -  Hk  +  Hk 
i/r  i  r  i  r 


where  T  is  the  observed  arrival  time  and  H  is  the  expected 
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(Herrin  1966)  travel  time  from  the  hypocenter  of  the  kth  event 
including  correction  for  ellipticities  but  not  for  station 
elevations . 

12.  A  fixed-point  zero  anomaly  indicates  that  no  readina 

was  made  at  the  subarray  for  that  event. 

13.  The  average  anomaly  at  subarray  i  of  N  recorded  events? 


/N 


for  the  defined  region. 

14.  Standard  deviation,  or  error  of  estimate,  at  the  ith 


subarray  for  N  observations: 


0  -J 

r 

,  _  ri 

Ji 

/ (N-l)  - 

1  1 

i 

-< 

■>»  1 

for  the  defined 

region. 

>5 


15.  Number  of  observations,  N,  at  station  i  for  the  defined 


region. 

16.  Total  number  of  epicenters  included  in  the  defined 


region. 

17.  Epicenter  latitude,  degrees  (USC&GS) ;  plus  north, 
minus  south. 

18.  Epicenter  longitude,  degrees  (USC&GS) ?  plus  east, 
minus  west. 

19.  Event  depth,  kilometers  (USC&GS) . 

20.  Event  origin  time,  hours,  minutes,  seconds  (USC&GS). 

21.  Standard  deviation,  or  error  of  estimate,  of  the  kth 
event  in  the  defined  region; 
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35 


/(L-l)> 


where  L  is  the  number  of  subarrays  recording  the  kth  event  not 
including  the  reference  subarray  R. 

22.  Average  error,  or  bias,  of  the  kth  event; 


L 

E.  =  V 


...  (Ak,  -  A.  )/L 
1=1  l/r  i/r 


where  L  is  the  number  of  subarrays  recording  the  kth  event  not 
including  the  reference  subarray  R. 

23.  Number  of  subarrays,  L,  recording  the  kth  event,  not 
including  the  reference  subarray  R. 

The  results  in  Report  204  are  arranged  first  by  general 
direction  (beginning  with  the  northwest  and  going  clockwise  and 
second  by  increasing  epicentral  distance  within  each  directional 


group  (See  Table  1) . 

The  following  observations  similar  to  those  in  a  previous 
report  using  the  Jeffreys-Bullen  travel-time  tables  (Chiburis, 
1966) ,  are  made  concerning  the  LASA  travel-time  anomalies  using 
the  Herrin,  November  1966  tables: 

1.  The  anomaly  variations  between  regions  measured  at  a 

single  subarray  are  as  high  as  1.66  seconds.  For  example,  sub¬ 
array  F4  has  an  average  anomaly  of  +0.74  sec  (N  =  19,  =  0.13) 

for  the  No.  Colombia  region  and  an  average  of  -0.92  sec  (N  =  12, 

a  .  =  0.07)  for  the  Solomon  Is.  region. 

F4 

2.  Subarrays  which  are  quite  near  the  reference  subarray 
A0  can  have  unusually  large  anomalies.  For  example,  subarray  D4 
is  located  30.75  km  from  A0  (center  instrument  to  center  instru¬ 
ment)  and  has  an  anomaly  relative  to  A0  of  +0.83  sec.  (N  =  7, 

a  =  0.06)  for  the  Dominican  Republic-Mona  Passage  region. 

D4 
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TABLE  I 


Distance-azimuth  ranges  by  region 


Direction 


No. 

cf  Distance  Azimuth 

Rgcri -ins  Range,  km  Range,  deg 


No. 

of 

Events 


Northwest 

North 

Northeast 

East 

East-southeast 

Southeast 

South 

Southwest 

West 

Undefined 

Continental  U.S 
Eastern  Is.  and 
Pacific  Ocean 
Miscellaneous 


15  2800-10900 

6  5300-11000 

9  4/00-10700 

4  5700-  9900 

3  4500-  5800 

13  3100-10800 

8  2200-10100 

2  9500-10900 

2  10600-11100 

3 

683-  2900 

4400-10800 

2100-9600 


292-329 

214 

340-016 

36 

018-071 

65 

075-103 

16 

112-121 

23 

132-170 

133 

158-192 

49 

238-248 

43 

259-274 

20 

101-281 

7 

184-253 

6 

334-007 

14 

3.  The  center  instrument  at  subarray  B2  is  only  7.50  km 

from  the  center  at  AO,  but  it  has  an  anomaly  of  -0.30  sec  (N  =  13, 

a  =  0.04)  for  the  Yugoslavia-Albania-Greece-Mediterranean  Sea 
B2 

Region.  This  result  suggests  that  the  time  anomalies  within  one 
subarray  (7  km  diameter)  are  far  from  negligible.  Signals  re¬ 
corded  within  a  subarray  can  be  significantly  misaligned  with 
anomalies  as  large  as  0.30  sec. 

4.  The  anomalies  are  not  slowly  varying  functions  of  either 

distance  or  azimuth.  For  example,  subarray  Fl  has  an  anomaly  of 

+0.22  sec  (N  =  7,  a  =  0,09)  computed  from  events  approaching 

Fl 

from  an  east-southeastly  direction  at  5100  km  distance  (Virgin- 

Leeward  Is.  region),  but  at  4600  km  (Dominican  Republic-Mona 

Subarray  F2,  on  the  other  hand,  has  an  anomaly  of  +0.12  sec 

(N  =  6,  =  0.08)  for  events  bearing  145°  at  4800  km  (So. 

,  o 

Central  America) ,  whereas  for  events  bearing  113  at  a  distance 

of  5100  km  (Virgin-Leeward  Is.  region)  the  anomaly  is  -0.57  sec 

(N  =  6,  o  =  0.06).  Hence,  the  anomaly  at  Fl  changes  by  0.62 
F2 

sec  in  a  distance  range  of  500  km,  and  at  F2  it  changes  by  0.69 
sec  in  an  azimuth  range  of  32°. 

5.  The  maximum  anomaly  range  observed  at  LASA  is  1.94  sec; 
average  anomaly  at  subarray  F2  is  -1.01  sec  for  events  occurring 
in  Rumania;  average  anomaly  at  subarray  El  is  +0.93  for  events 
from  the  No.  Colombia  region. 

6.  The  maximum  anomaly  range  for  one  particular  region 
(North  Atlantic  Ridge)  is  1.43  sec,  where  the  D4  anomaly  is 
+0.47  sec  and  the  F2  anomaly  is  -0.97. 

\ 
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Chiburis,  E.F. ,  1966,  "LASA  Travel-Time  Anomalies  for  Various 
Epicentral  Regions",  Seismic  Data  Laboratory  Report  No.  159, 
13  September. 

B.  Lateral  Variations  in  Crustal  Structure  Beneath  the 

Montana  LASA 

The  analysis  of  a  variety  of  geologic  and  geophysical 
data  indicates  that  significant  lateral  variations  in  crustal 
structure  exist  across  LASA  This  structural  complexity  is  in¬ 
ferred  from  observation  of  Rayleigh  wave  dispersion,  body-wave 
spectral  ratios,  P-wave  amplitude  and  travel-time  anomalies, 
seismic  refraction  profiles  and  gravity  and  magnetic  anomalies. 

To  explain  these  results,  LASA  is  divided  into  an  eastern  and  a 
western  sector.  An  attempt  is  made  to  explain  the  observations 
in  each  sector  in  terms  of  crustal  models  consisting  of  3  to  5 
distinct  layers.  The  observations  show  that  the  maximum  lateral 
change  in  structure  across  LASA  takes  place  in  a  NE-SW  direction. 

The  most  striking  feature  of  the  evidence  presented  is 
its  consistency.  The  repeated  appearance  of  the  NW-SE  trend 
suggests  that  this  is  fundamental  to  the  structure  beneath  the 
array.  Only  in  the  distribution  of  travel-time  anomalies  is  it 
not  readily  discernible.  In  the  following  discussion  we  shall 
attempt  to  explain  our  observations  in  terms  of  a  model  which 
may  form  the  basis  for  future  investigations. 

From  the  Rayleigh  wave  dispersion  data,  supplemented  by 
the  evidence  from  the  spectral  ratios  and  the  aeromagnetic  map, 
it  is  evident  that  as  a  starting  point  LASA  may  be  divided  into 
two  sectors  by  the  line  A-A'  (See  Figure  2).  To  the  east  of 


Figure  2.  Location  of  the  Montana  LASA 


of  this  line  the  observed  phase  velocities  are  higher  than  those 
to  the  west.  Of  the  crustal  models  considered,  U.W. 3  best  fits 
the  observed  data  for  the  eastern  sector,  whereas  USGS3  best 
fits  the  observed  data  for  the  western  sector  of  LASA.  Some 
important  implications  of  this  division  are: 

1.  Seismic  velocities  in  the  eastern  sector  are  sig¬ 
nificantly  higher  than  those  in  the  western  sector; 

2.  The  total  thickness  of  the  crust  is  greater  in 
the  east  than  in  the  west; 

3.  The  models  indicate  a  reduction  in  the  number  of 
major  layers  within  the  basement,  from  three  in  the  east  to  two 
in  the  west. 

The  observed  decrease  in  phase  velocities  to  the  west 
of  A-A'  is  in  apparent  contradiction  to  the  overall  reduction 
in  thickness  of  the  crust  inferred  from  the  seismic  refraction 
results.  To  explain  this  it  is  necessary  to  postulate  the 
presence  either  of  dipping  layers  within  the  crust,  or  lateral 
changes  in  velocity  within  particular  layers,  or  both.  Dipping 
layers  have  been  identified  in  the  seismic  refraction  studies 
(See  Figures  3  and  4) ;  furthermore,  the  implied  reduction  in 
the  number  of  layers  within  the  basement  may  be  explained  in 
this  manner. 

The  nature  of  the  boundary  between  the  two  sectors  is 
problematical.  The  geologic  data  do  not  indicate  the  presence 
of  large  fault  or  fault  zone  within  the  sedimentary  sequence 
(i.e.,  the  uppermost  3  km)  of  the  crust.  This  is  emphasized 
by  the  dispersion  results,  which  indicate  that  the  amount  of 
lateral  refraction  which  Rayleigh  waves  undergo  on  crossing  LASA 


-7- 


Figure  4.  The  USGS  Crustal  Model.  (After  Borcherdt  and  Roller,  1967) 


increases  with  depth,  at  least  to  about  50  km.  This  implies 
that  the  boundary  lies  at  some  depth  within  the  crust,  as  has 
been  postulated  by  Zeitz  (personal  communication,  1967)  on  the 
basis  of  the  magnetic  data. 

Moreover,  it  is  possible  that  the  boundary  is  gradational 
rather  than  abrupt.  The  dispersion  results  indicate  only  that 
a  change  takes  place.  We  must  await  the  detailed  interpretation 
of  the  gravity  and  magnetic  maps  before  speculating  on  the  nature 
of  the  boundary  itself. 

So  far  we  have  sought  only  to  explain  the  dispersion  results. 
The  travel-time  anomalies  must  also  be  explained  in  terms  of  a 
crustal  model.  Of  the  explanations  put  forward  to  date,  three 
will  be  discussed  here;  bearing  in  mind  that  any  explanation 
must  satisfactorily  account  for  the  "synclinal"  feature  reported 
by  Sheppard  (1967)  and  substantiated  by  our  results. 

To  explain  the  synclinal  feature  Sheppard  (1967)  has  pro¬ 
posed  two  alternate  hypotheses.  They  are;  . 

1.  A  relative  thickening  of  the  earth's  crust  by  soms 
10  km  under  subarray  B4  with  respect  to  F2  and  F4-  This  would 
correspond  to  a  five-degree  to  six-degree  slope  in  the  Moho 
from  F4  to  F2  and  from  F2  to  B4. 

2.  A  thickening  of  the  3.0  km/sec  (sedimentary)  layer 
by  3  km  under  subarray  B4. 

Neither  of  these  hypotheses  are  supported  by  our  results. 

If  the  Moho  dips  by  as  much  as  5  degrees  then  the  division  of 
the  array  into  northern  and  southern  sectors  should  yield  rela¬ 
tively  lower  phase  velocities  in  the  north  compared  to  the  south, 
for  the  Greenland  Sea  event.  As  was  indicated  'previously,  the 
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phase  velocities  for  the  southern  sector  were  higher  than  those 
for  the  northern  sector.  Also  the  division  of  LASA  into  an 
inner  and  outer  group  of  subarrays  indicated  no  strong  thickening 
under  the  center  of  the  array.  From  a  consideration  of  the  well 
log  data  and  from  the  geologic  sections  presented  by  Brown  and 
Poort  (1965) ,  there  is  no  evidence  to  suggest  an  increase  of  3  km 
(nearly  100%)  in  the  sedimentary  sequence  anywhere  under  LASA. 

The  third  hypothesis,  proposed  by  Fairborn  (1966) ,  suggests 
that  the  azimuthal  dependency  of  the  travel-time  anomalies  re¬ 
quires  horizontal  velocity  gradients  in  the  crust  and  upper  mantle. 
Our  data  and  Sheppard's  (1967)  indicates  that  horizontal  velocity 
gradients  are  insufficient  in  themselves  to  account  for  the  ob¬ 
served  anomalies.  Both  the  seismic  refraction  data  and  the  Rayleigh 
wave  dispersion  data  indicate  that  the  average  velocity  of  a  55  km- 
thick  section  under  LASA  decreases  from  east  to  west.  Such  a  de¬ 
crease  in  velocity  would  result  in  systematically  later  arrivals 
towards  the  west  of  the  array  if  horizontal  velocity  gradients 
alone  are  responsible  for  the  observed  travel  time  anomalies.  From 
our  observations  and  Sheppard's,  the  later  arrivals  define  a 
"synclinal"  feature  to  the  NE  of  AO.  This  suggests  that  two  of 
the  more  important  factors  bearing  on  the  distribution  of  the  ob¬ 
served  travel  time  anomalies  are: 

1.  dipping  layers  within  the  crust? 

2.  the  pre-Cretaceous  structure  of  the  Paleozoic  rocks 
and  possibly  the  upper  part  of  the  basement  also. 
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Moreover,  since  large  amplitudes  are  consistently  recorded 
at  the  center  of  LASA,  it  seems  that  the  anticlinal  structures 
which  surround  the  array  giva  rise  to  a  lens  effect  which  re¬ 
sults  in  the  focusing  of  body  waves  toward  the  center  of  the 
array.  However,  more  observations  are  needed  to  substantiate 
that  this  occurs. 

From  the  analysis  of  the  geologic  and  geophysical  data 
we  conclude: 

1.  Lateral  variations  exist  in  the  structure  of  the 
earth's  crust  beneath  LASA.  The  geologic,  aeromagnetic  and 
surface  wave  dispersion  data  show  that  this  variation  is  within 
the  basement  at  a  depth  greater  than  three  kilometers.  The 
maximum  change  takes  place  in  a  NE-SW  direction  (i.e.,  perpen¬ 
dicular  to  the  dominant  NW-SE  structural  trend) . 

2.  The  observations  are  best  explained  by  dividing 
LASA  into  an  eastern  and  a  western  sector  by  a  line  trending 
N20°W  through  Miles  City,  Montana.  For  the  eastern  sector,  the 
average  seismic  velocities  of  the  crust  are  higher,  the  crust  is 
thicker,  and  Model  U.W. 3  best  accounts  for  the  observed  disper¬ 
se'  n.  For  the  western  sector,  the  velocities  are  lower,  the 
crust  is  thinner,  and  model  USGS3  best  accounts  for  the  observed 
dispersion  results. 

3.  The  individual  layers  within  the  crust  are  dipping, 
though  not  in  the  same  direction. 

To  refine  the  model  further  it  is  recommended  that: 

1.  The  individual  travel-time  anomalies  at  each  station 
be  obtained  as  a  function  of  distance  and  azimuth; 


-10- 


2.  More  events  be  examined.  In  particular,  more  spectral 
ratios  must  be  obtained  to  determine  an  average  structure  under 
each  subarray  and  then  Rayleigh  waves  from  events  at  more  azimuths 
and  distances  should  be  analyzed  to  establish  better  the  attitude 
and  thickness  of  the  intermediate  crustal  layers. 
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C.  High-Resolution  Frequency-Wavenumber  Spectra 

Burg  (1967)  has  shown  that  high-resolution  power  spectra 
can  be  computed  from  relatively  short  correlation  functions  by 
first  least-squares  predicting  larger  correlation  lags  from  those 
already  known.  The  resulting  high-resolution  power  spectrum  is 
inversely  proportional  to  the  power  spectrum  of  the  corresponding 
prediction  error  operator.  However,  with  the  advent  of  the 
Cooley-Tukey  fast  Fourier  transform  algorithm  and  related  methods 
for  estimating  power  spectra,  the  high-resolution  technique  for 
ordinary  power  spectra  is  no  longer  economical.  It  is  now  feasible 
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to  compute  power  spectra  as  if  they  were  Fourier  transforms  of 
the  100%-lag  correlation  functions;  various  kinds  of  smoothing 
can  then  be  applied  to  reduce  the  spectra  to  the  desired  length 
and  stability. 

We  feel  that  the  most  promising  application  of  the 
high-resolution  procedure  lies  in  estimating  frequency-wavenumber 
spectra „  Here  the  length  of  the  spatial  correlation  functions, 
which  ultimately  determines  the  resolution,  is  limited  by  the 
size  of  the  array.  The  high  resolution  method  is  therefore  the 
appropriate  way  to  obtain  frequency -wavenumber  spectra  which  are 
sharp  in  wavenumber  space  from  data  recorded  over  small  arrays. 

It  should  be  noted  that  computing  high-resolution  fre¬ 
quency-wavenumber  speccra  wi th  our  method  is  only  negligibly  more 
expensive  than  computing  ordinary  frequency-wavenumber  spectra. 
Both  techniques  require  computing  the  spectral  matrix.  The  only 
other  computations  required  by  the  high-resolution  technique  are 
the  addition  of  some  white  noise,  one  Hermitian  matrix  inversion, 

and  one  Hermitian  matrix  multiplication. 

Two  events  and  samples  of  noise  immediately  preceding 
them  were  selected  as  examples  of  high-resolution  frequency- 
wavenumber  data  processing.  The  first  event  and  its  corresponding 
noise  sample  were  used  to  illustrate  the  application  of  the  tech¬ 
nique  to  a  two-dimensional  surface  array,  in  this  case  the  WMSO 
array,  while  the  remaining  data  were  used  to  illustrate  results 
from  a  one-dimensional  vertical  array,  the  APOK  deep  well. 

Ordinary  frequency-wavenumber  spectra  of  both  pairs 
of  noise  samples  and  signals  are  shown  in  Figures  5,  7,  9a  and 
10a.  Figures  5  and  9a  are  the  noise  samples  and  Figures  7  and 
10a  are  the  signals.  The  WMSO  noise  seems  to  be  coming  from  the 
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northeast  with  a  velocity  ranging  between  3  and  6  km/sec-  This 
characteristic  behavior  has  been  noted  by  other  investigators 
(Rizvi  et  a_l .  ,  1967).  On  the  other  hand,  the  WMSO  signal 
appears  to  be  coming  from  the  south  or  southeast  with  a  velocity 
somewhere  between  16  and  20  km/sec.  The  APOK  noise  sample  dis¬ 
plays  a  tilted  pattern  also  characteristic  of  the  site.  One 
explanation  for  this  phenomenon  is  suggested  by  the  geological 
structure  in  the  vicinity  of  the  well  (Geotechnical  Corp.,  1964), 
which  consists  of  Paleozoic  intrusives  and  sediments  dipping 
30-40°  towards  the  northeast.  The  APOK  signal  is  very  strange 
indeed,  composed  of  a  wave  travelling  down  the  well  with  an 
apparent  vertical  velocity  of  approximately  3.5  km/sec  and  a 
horizontally  travelling  wave  with  infinite  vertical  '  elocity. 

One  possible  ejgjlanation  of  this  result  is  again  suggested  by 
the  dip  of  the  structure  (Sax,  1967) .  The  incoming  P  wave  might 
be  refracted  by  the  tilted  interface  so  that  the  wave  front 
would  be  a  leaking  surface  wave  when  it  reached  the  well. 

High-resolution  frequency-wavenumber  spectra  are  shown  for 
all  these  data  samples  in  Figures  6,  8,  9b  and  10b.  These  were 
computed  with  a  signal-to-noise  ratio  of  2,  which  seemed  to  be 
a  good  balance  between  high  resolution  and  excessive  distortion. 
For  the  WMSO  examples,  the  technique  decreased  the  width  of  the 
-3db  contour  by  more  than  a  factor  of  two  for  both  the  cases  of 
noise  and  signal.  The  results  are  more  difficult  to  assess  for 
the  vertical  array  data  because  of  the  "whitening"  in  frequency 
done  by  the  filters.  However,  for  the  noise  sample  at  0.25  cps, 
a  decrease  by  a  factor  of  two  in  the  width  of  the  -3db  contour 
can  again  be  seen. 
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A  marked  disadvantage  and  a  difficulty  still  to  be  solved 
in  the  application  of  this  technique  is  illustrated  in  Figure  10b. 
The  high— re solution  technique,  when  applied  to  short  samples  of 
data,  can  produce  spurious  peaks  in  the  frequency-wavenumber 
spectra.  These  are  due  to  the  limited  smoothing  done  when  com¬ 
puting  spectra  from  short  samples  of  data.  Figure  10b  shows 
several  of  these  peaks.  This  is  regarded  as  an  unrecoverable 
distortion  and  in  a  certain  sense,  an  inevitable  result  of  using 
a  high-gain  procedure. 
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L.  Best  Linear  Unbiased  Estimation  for  Multivariate 

Stationary  Processes 

Many  problems  in  the  area  of  applied  time  series  analysi 
can  be  formulated  and  solved  by  extending  and  generalizing  well 
known  techniques  from  the  Isssical  theory  of  the  multivariate 
linear  hypothesis.  The  ft  in'ogies  between  some  of  the  physical 
models  expressed  in  terms  of  signals  propagating  across  arrays 
and  the  usual  analysis  of  variance  models  involving  various  kinds 
of  effects  are  striking,  particularly  when  it  can  be  assumed  that 
the  multivariate  stochastic  process  in  question  is  weakly  sta¬ 
tionary.  In  the  stationary  case,  simplicity  in  exposition  as 
well  as  economy  in  computation  result  from  performing  the 
analysis  in  the  frequency  domain  using  the  properties  of  the 
Fourier  transform. 

As  an  example  of  a  simple  estimation  problem  suppose 
that  a  multivariate  stochastic  process  fc\(t),  j=l,2,...,n, 

-  oo  <t  <<»}  consists  of  a  signal  which  is  identical  for  each  j 
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and  a  weakly  stationary  multivariate  noise  process.  More 
specifically  we  assume  that 

Y_,(t)  =  s(t)  +  n_,(t)  j=l,2,...,n 

—  co<  <oo 

In  this  model  s(t)  is  regarded  as  a  fixed  unknown  signal  with 
n.(t)  a  weakly  stationary  zero  mean  multivariate  noise  process 
with  a  correlation  function  given  by 

Rjk(t-t')  "  EVt)nk(t,)  = 


where  we  assume  that  0^(w)  is  cross  spectral  density 

matrix  of  the  noise  process.  The  notation  is  simplified  if  we 
use  continuous  parameter  processes  and  the  results  apply  equally 
well  to  the  discrete  parameter  sampled  data  version.  By  a  best 
linear  unbiased  estimate  (BLUE)  of  the  signal,  say  s(t),  is  meant 


the  unbiased  linearly  filtered  version  of  the  data  (Es(t)  =  s(t)) 

A  2 

which  has  the  smallest  variance  (E(s(t)  -  s(t))  =  min.).  This 

problem  is  analogous  to  estimating  the  mean  of  a  set  of  correlated 
random  variables  and  has  been  previously  considered  by  Kelly  and 
Levin  (5) ,  and  Capon,  Greenfield,  and  Kolker  (1) .  The  inherent 
advantage  of  this  model  over  the  Wiener  approach  (for  example  see 
Papoulis,  (8)  is  that  it  is  not  necessary  to  know  in  advance  the 
functional  form  the  signal  in  order  to  minimize  the  mean  square 
error  or  variance  of  the  signal  estimate. 

The  general  solution  for  BLUE  estimators  given  by  H(to)  = 
(X*  I  XX)  1X*  J  ^  (w)  requires  that  the  spectral  matrix  of  the 
noise  be  known;  although  this  is  seldom  the  case,  good  estimates 
for  l  can  sometimes  be  obtained  by  using  a  sample  of  noise  when 
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the  regression  coefficients  3 ^ ( t)  are  not  present.  In  the  next 
section  we  show  how  the  general  linear  model  through  an  appropriate 
choice  for  X_k(t)  can  be  specialized  to  a  number  of  interesting 
cases  occurring  in  practical  applications. 

The  theoretical  results  are  easily  extended  to  discrete 


time  parameter  processes  by  replacing  integrals  with  summations 
and  the  infinite  frequency  range  by  (-tKw<tt  )  but  in  the  applica¬ 
tion  to  finite  time  sampled  data  certain  approximations  are  neces¬ 


sary.  We  must  decide  first  whether  the  analysis  proceeds  more 
reasonably  and  economically  in  the  time  domain  or  the  frequency 
domain.  We  have  chosen  the  frequency  domain  for  several  reasons. 
First,  the  restrictions  imposed  by  certain  signal  models  must 
allow  for  time  delays  as  long  as  300  digital  points.  This  means 
that  the  time  domain  analogue  of  the  matrix  product  (X*  £  X) 
becomes  extremely  large  if  a  512  point  time  delay  preserving 
filter  is  required.  Hence,  the  frequency  domain  approach  be¬ 
comes  much  faster  as  the  analysis  can  be  performed  separately 
at  each  frequency.  In  addition,  it  has  been  noted  in  Capon, 
Greenfield,  and  Kolker,  (1)  that  short  filters  designed  in  the 
time  domain  are  more  sensitive  to  slight  stationarity  and  signal 
model  perturbations.  Reference  1  also  contains  several  examples 
illustrating  time  and  frequency  domain  computations  and  concludes 
that  the  loss  due  to  using  a  two  sided  infinite  lag  filter  rather 
than  the  physically  realizable  filter  obtained  in  the  time  domain 
is  small. 

The  approximations  used  in  the  examples  below  are  based 
on  the  use  of  the  finite  Fourier  transform  and  involve  replacing 
the  data  and  filter  by  aliased  versions  of  the  optimum  infinite 
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two-sided  Fourier  transforms.  References  2  and  7  contain  excellent 
material  describing  the  theory  and  application  of  the  finite 
Fourier  transform.  The  effects  of  aliasing  are  reduced  by  choosing 
an  appropriate  sampling  interval  in  time  based  on  the  highest 
frequency  observable  in  the  data.  In  the  two  examples  given  below 
the  signals  are  limited  to  the  band  from  0  to  10  cycles  per  second 
(cps)  with  the  major  frequency  content  less  than  5  cps.  Hence, 
the  data  can  be  sampled  with  A  t  =  .05  seconds  yielding  an 
aliasing  frequency  of  10  cps.  The  data  samples  are  long  with  the 
first  being  150  seconds  (3,000  pts)  and  the  second  60  seconds 
(1200  points)  so  that  -<=°<  t<°°  is  approximately  valid.  Finally, 
the  accurticy  of  the  various  approximations  employed  can  be  eval¬ 
uated  from  simulated  examples  similar  to  those  given  below. 

Example  1 

Here  we  consider  the  model  given  by 

Y .  ( t)  =  s.  (t)  +  s  (t-T  )  +  n  (t) 

3  J-  *  j  j 

j  =  1 , 2 ,  .  . .  ,  20 

This  model  is  equivalent  to  two  interfering  signals  recorded 
at  the  same  time. 

Fig  are  11  shows  the  first  five  channels  of  data  con¬ 
structed  to  conform  with  the  above  equation.  Two  seismic  signals 
and  smoothed  white  noise  uncorrelated  from  channel  to  channel  were 
added  in  the  same  proportions  to  produce  twenty  channels.  Five 
of  the  twenty  channels  are  displayed  in  Figure  11.  The  general 
signal  estimation  procedure  was  to  compute  the  matrix  product 
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Yj [♦)  Y^t)  Y^t)  Y4(t)  Y5(t;i 


Figure  11  Separation  of  Two  Mixed  Signals 
By  Minimum  Variance  Filters 


(X*X)-1X*(w  )  at  each  of  256  frequencies.  The  singularity  at 

w=  0  was  eliminated  by  tapering  the  frequency  response  functions 

down  to  zero  at  0  cps.  A  fast  Fourier  transform  subroutine 

(see  McCowan,  7)  applied  to  the  frequency  response  matrix  H(w) 

produced  the  impulse  response  functions  shown  in  Figure  12. 

Note  that  the  net  results  of  the  filters  are  to  reinforce  and 

sum  the  aligned  signal  s^t)  while  cancelling  the  second  signal 

s  (t)  which  appears  at  the  given  time  delays.  The  reverse  holds 

when  estimating  s2(t)  with  the  unaligned  signal  reinforced  and 

the  aligned  signal  cancelled.  The  signal  estimates  are  formed 

using  a  high  speed  convolution  subroutine  which  calculates  a 

matrix  product  in  the  frequency  domain  and  then  re-transforms 

the  result  to  obtain  time  domain  estimates  for  the  regression 

functions.  Note  that  in  the  application  of  the  convolution  each 

Y  (t)  is  available  only  for  the  finite  interval  T  and  we  assume 
j 

that  Y .  ( t)  =0  for  t  &  T.  The  true  signals  s^t)  and  s2(t)  are 
displayed  along  with  the  estimates  s^t)  and  s2(t)  in  Figure  11. 
The  generation  of  the  filters  at  256  frequencies  or  512  time 
points  for  20  channels  took  eight  minutes  of  CDC  1604  time  with 
the  subsequent  convolution  requiring  thirty  minutes.  However, 
by  combining  the  operations  in  a  different  order  the  total 
running  time  has  been  reduced  to  eight  minutes. 

Example  2 

As  an  example  of  a  more  complicated  model  consider  the 
multivariate  process  recorded  on  a  vertical  array  shown  in 
Figure  13.  Each  level  or  channel  conforms  artificially  to  the 

model, 
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Y.(t)  =  6.  (t)  +  t*  r00  x.  (t-u)  0  (u)  du  +  n  (t) 

J  1  J  -oo  m  3 


(j  =  1 . 5) 


The  data  representing  a  teleseismic  signal  and  several 

Rayleigh  modes,  was  generated  by  adding  noise  and  the  known 

functions  g.(t)  in  Figure  13  to  the  data  using  the  transfer 

functions  X.,  (t)  .  The  real  valued  frequency  responses  (X  (f)m 
jk  JK 

f  =  w  /2  cps)  are  shown  in  Figure  14  for  the  range  0  to  2  cps. 
The  response  functions  were  low  pass  filtered  with  a  cutoff  at 
2  cps:  hence,  there  is  no  significant  frequency  response  above 
that  value.  The  optimum  filters  in  this  case  were  200  points  or 
10  seconds  long  and  were  computed  by  evaluating  the  matrix 
product  (X*X)-1X*(w  )  at  200  frequencies.  (The  fast  Fourier 
transform  was  not  used  here  so  the  number  of  data  points  is  not 
a  power  of  two) .  The  resulting  matrix  of  filter  coefficients 
was  convolved  with  a  multivariate  process  to  generate  the 
estimates  shown  in  Figure  13.  Comparison  of  these  estimated 
regression  functions  with  the  true  regression  functions  in¬ 
dicates  that  the  procedure  again  produces  reasonable  estimates. 
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E.  Design  and  Evaluation  of  Certain  Multichannel  Filters 
The  first  phase  of  the  SDL  multichannel  filter  project, 
namely,  the  design  and  evaluation  of  certain  multichannel  filters 
has  been  concluded.  Experimental  results  obtained  by  applying 
the  various  multichannel  filters  to  actual  seismic  array  data 
are  used  to  reach  certain  conclusions  about  the  performance  of 
the  various  multichannel  filters. 

The  practical  computation  of  many  of  these  filter  opera¬ 
tors  relies  heavily  on  the  fast  Fourier  transform  algorithm  of 
Cooley  and  Tukey,  without  which  many  of  the  calculations  would 
not  even  be  possible.  An  earlier  report  (McCowan,  1966)  dealt 
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with  the  transform  method  and  its  associated  techniques  and  the 
present  report  can  be  regarded  as  a  logical  development  of  the 
ideas  presented  there. 

The  standard  of  comparison  used  to  grade  multichannel 
filter  performance  is  the  unbeamsteered  sum  of  the  data  channels. 
This  was  chosen  instead  of  the  beamsteered  or  phased  sum  simply 
as  a  matter  of  convenience.  Straight  summation  is  itself  an 
optimum  processor  when  the  signals  are  aligned  and  the  noise  is 
uncorrelated.  Any  departures  from  this  model  either  in  signal 
misalignment  of  noise  correlation,  will,  of  course,  detract  from 
the  theoretically  best  performance,  which  is  the  square  root  of 
the  number  of  sensors  increase  in  signal  to  noise  (S/N)  ratio 
improvement.  This  is  where  the  multichannel  filter  is  applicable; 
the  design  equations  are  a  theoretical  statement  of  these  depart¬ 
ures  in  signal  and  noise  models.  Filters  satisfying  these  equa¬ 
tions  should  provide  an  increase  in  S/N  ratio  improvement  over 
the  otherwise  optimum  processor,  the  straight  sum. 

The  data  used  in  the  experimental  part  of  this  report 
were  from  an  Aleutian  earthquake  recorded  by  ten  sensors  in  the 
Dl  subarray  of  the  Montana  LASA.  Ten  sensors  were  used  instead 
of  the  original  twenty-five  available  in  the  Dl  subarray  to 
reduce  the  computation  time  and  to  simplify  the  procedure  of 
tabulating  results.  This  particular  set  was  selected  to  provide 
a  reasonably  geometric  array  with  good  element  separation.  Cali¬ 
brations  were  computed  in  the  us-ual  manner  and  all  traces  were 

corrected  to  true  ground  motion. 

For  some  parts  of  the  study,  the  traces  were  prefiltered 
with  the  "SDL"  bandpass  filter.  This  filter  is  a  phaseless  four- 
pole  Butte.rw’orth  recursive  filter.  It  was  chosen  to  provide  a 
balance  between  excessive  signal  distortion  at  one  end  and  good 
noise  reduction  at  the  other. 


/ 


As  in  previous  SDL  da ta-processing  reports,  the  defini¬ 
tion  of  signal  amplitude  will  be  one-half  the  maximum  peak-to- 
trough  amplitude  occurring  in  the  first  three  cycles  of  the  P 
phase.  Noise  is  defined  as  the  root-mean-square  amplitude 
measured  over  some  specified  interval.  The  S/N  ratio  is  the 
quotient  of  these  two  quantities.  Frequently  results  will  be 
expressed  in  db  for  which  the  usual  expression  will  be  used, 

20  log  (RATIO) .  When  the  quantities  to  be  expressed  in  terms 
of  db  have  the  dimensions  of  power,  this  definition  has  been 
changed  to  10  log. _ (RATIO)  in  order  to  preserve  the  proper  am- 
plitude  ratio. 

The  results  of  this  study  are  presented  in  Table  2 
which  lists  the  actual  improvement  in  S/N  ratio  over  the  un- 
prefiltered  straight  sum  for  the  various  multichannel  filters. 
Simply  bandpass  filtering  the  sum  gets  15  db  improvement  in  S/N 
ratio  at  practically  no  expense  of  computer  time,  so  the  real 
comparison  is  with  this  reference  level.  The  amplitude  of  the 
noise  was  computed  for  the  measured  noise  multichannel  filter 
outputs  over  an  interval  of  200  seconds,  including  the  fitting 
interval  and  100  seconds  outside  the  fitting  interval.  The 
maximum-likelihood  noise  amplitude  was  also  computed  over  an 
interval  of  200  seconds,  but  this  time  including  the  fitting 
interval  and  50  seconds  outside  it.  Therefore,  these  noise 
levels  represent  an  average  of  the  rms  noise  amplitude  inside 
fitting  interval  and  outside  it.  This  problem  does  not,  of 
course,  arise  with  the  straight  sums  and  theoretical  isotropic 
processor  outputs. 
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TABLE  2 


MULTICHANNEL  FILTER  REPORT 


S/N  Ratio  Improvement 

Filter 

Data 

Signal  Model 

db  Improvement 

Unphased  Sum 

UBPF 

— 

0 

Theoretical  MCF 

UBPF 

10  km/ sec 

7 

Theoretical  MCF 

UBPF 

Infinite 

4.8 

Maximum  Likeli¬ 
hood 

UBPF 

— 

15 

Measured-Noise 

MCF 

UBPF 

10  km/sec  . 

3 

Measured-Noise 

MCF 

UBPF 

Infinite 

-6 

Unphased  Sum 

BPF 

— 

15 

Theoretical  MCF 

BPF 

10  km/sec 

16.4 

Theoretical  MCF 

BPF 

Infinite 

16.1 

Maximum  Likeli¬ 
hood 

BPF 

— 

17 

Measured-Noise 

MCF 

BPF 

10  km/sec 

10.3 

Measured-Noise 

MCF 

BPF 

Infinite 

5 

The 


We  can  make  several  comparisons  in  this  table, 
first  concerns  the  amount  of  frequency  filtering  that  is  done. 
The  illusory  gains  from  the  maximum-likelihood  filter  operating 
on  un~pref iltered  noise  come  about  because  the  least-squares 
equation  does  not  constrain  the  filter  to  do  otherwise.  Many 
degrees  of  freedom  in  the  filter  weights  are  misused  doing 
what  could  be  done  with  just  elementary  processing,  i.e.,  band¬ 
pass  filtering.  Apparently  all  but  2  db  of  the  original  15  db 
can  be  accounted  for  by  bandpass  filtering.  On  the  other  hand, 
the  design  equation  for  the  measured-noise  multichannel  filter 
specifies  as  a  desired  signal  spectrum  one  of  the  noise  spectra. 
Therefore,  some  frequency  filtering  is  expected.  The  design 
equation  for  the  theoretical  multichannel  filter  assumes  a  flat 
spectrum  for  the  desired  signal  output.  However,  the  filter 
which  satisfies  this  equation  is  only  the  optimum  filter  in  the 
least-squares  sense,  so  we  should  expect  some,  but  not  much, 
frequency  filtering.  This  is  demonstrated  in  the  table.  In 
summary,  the  maximum  likelihood  multichannel  filter  is  the 
worst  regarding  frequency  filtering  and  the  theoretical  multi¬ 
channel  filter  is  the  best. 

Secondly,  in  both  cases,  the  theoretical  and  the 
measured  noise  filters  designed  with  a  disc  signal  model 
(V  >_  10  km/sec)  performed  better  than  the  infinite-velocity 
filters.  This  fact  is  true  in  general  and  is  a  direct  result 
of  not  beamsteering  the  array.  However,  other  results  dis¬ 
cussed  below  suggest  that  this  gain  occurs  only  in  a  frequency 
band  higher  than  the  signal  band.  For  the  prefiltered  theoret¬ 
ical  multichannel  filter,  designed  on  noise  where  this  energy 
has  been  removed,  this  difference  only  amounts  to  0.3  db. 
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Finally,  as  regards  isotropic  processor  multichannel 
filters,  the  theoretical  processor  is  always  better  than  the 
measured-noise  processor.  This  is  a  strange  and  completely 
unexpected  result,  probably  arising  from  the  fact  that  the 
noise  seems  to  overlap  the  signal  model  in  k  space.  For  a 
measured-noise  multichannel  filter  this  represents  a  contra¬ 
diction  of  specification  in  the  least-squares  design  equations 
which  could  lead  again  to  a  misuse  of  the  degrcses  of  freedom 
available  in  the  filter  weights.  The  theoretical  multichannel 
filter  is  not  computed  from  actual  data  and  hence  is  not  af¬ 
fected. 

The  conclusions  listed  below  are  based  upon  results 
obtained  from  processing  only  one  event.  In  view  of  this 
they  should  be  accepted  with  due  caution. 

1.  Designing  multichannel  filters  on  un-pref iltered 
noise  gives  an  illusory  S/N  ratio  improvement  which  is  not 
attained  with  similar  multichannel  filters  designed  or:  pre¬ 
filtered  noise. 

2.  Isotropic  processors  can  be  effective]. y  applied 
to  unbeamsteered  arrays,  provided  a  finite  velocity  signal 
model  is  used. 

3.  Measured-noise  isotropic  processors  are  not  ef¬ 
fective  in  improving  S/N  ratio  when  the  noise  has  large  high- 
velocity  components.  In  this  case  the  theoretical  multichannel 
filter  is  the  superior  processor. 

4.  With  the  exception  of  the  maximum-likelihood 
filter  evaluated  for  performance  within  its  fitting  interval, 
only  the  theoretical  multichannel  filter  showed  a  gain  over  the 
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prefiltered  straight  sum.  That  gain  was  slightly  more  than 
1  db,  which  is  not  regarded  as  significant. 

5.  One  drawback  of  the  maximum-likelihood  filter  is 
its  unnecessarily  complex  response  in  k-space.  It  does  not 
appear  that  modifying  the  constraint  will  cure  this. 
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F.  Digital  Computer  Programs  for  the  Design  and  Evaluation 
of  Multichannel  Filters 

The  following  is  a  list  of  program  write-ups  for  the 
multichannel  filter  pro<^  "am  set  written  under  the  multichannel 
filter  project  at  the  SDL.  All  of  these  programs  were  verified 
using  write-ups  which  are  now  as  free  from  errors  as  possible. 
This  set  of  programs  should  provide  a  reasonably  complete 
capability  to  analyze  signals  and  noise,  design  multichannel 
filters  to  enhance  signals  and  suppress  noise,  evaluate  the 
performance  of  these  filters,  and  prepare  the  punched  paper 
tape  which  inputs  the  multichannel  filters  into  the  Texas 
Instruments  Digital  converter. 
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An  alphabetical  list  of  the  programs  described  in 
SDL  Report  210  along  with  corresponding  COOP  identification 


numbers  is  shown  below: 


G621 

COLYTUKY 

M2  21 

MCFONPT 

G627 

COHERENCY 

G630 

MULTICOH 

G638 

CONFIL 

G631 

PARTLCOH 

G624 

FKSPTRUM 

G636 

PLOTFK 

G632 

HEFALUMP 

G628 

PREDICT 

G640 

ISOFIL 

G635 

SPECAVE 

G639 

LSTCHNCE 

G625 

TWX 

Z124 

MAXLIK 

G629 

VFKSPTRUM 

M2  20 

MCFONPT 

Response  of  Several 

Vertical 

Array  Processors 

The  SDL  frequency-wavenumber  analysis  program  VFKSITRM 
(1966)  was  used  to  calculate  the  response  of  several  simple 
processors  for  vertical  arrays.  The  impulse  responses  of  the 
filters  as  applied  to  each  channel  were  subjected  to  frequency- 
wavenumber  analysis.  This  implied  an  equal  and  simultaneous 
impulsive  input  to  each  channel. 

Each  processor  involves  a  beamed  sum  of  P-waves  based 
on  normal  moveout  times  for  vertical  compressional  waves.  The 
multichannel  deghoster  and  fan  filter  are  followed  by  other 
operations  designed  to  remove  signal  distortion  or  improve  the 
<3ead-band  response  of  the  processor.  The  multichannel  deghoster 
estimates  the  signal  waveform  at  the  surface  without  regard  to 
correlated  background  noise.  The  fan  filter  (1967)  ideally 
passes  waves  between  a  specified  pass-band  of  limiting  vertical 
phase  velocities  and  rejects  energy  outside  of  the  band. 
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We  examined  several  vertical  array  geometries.  These 
are  specified  by  uphole  times  based  on  the  arrival  of  the  up- 
going  vertical  compressional  wave  at  each  instrument  in  the 
array  subtracted  from  the  arrival  time  at  the  surface.  A 
vertical  phase  velocity  of  6  km/sec  was  used  for  the  f-k  analysis. 

Table  3  correlates  the  array  geometry,  type  of  processor 
and  figure  number  from  the  original  report.  Table  4  illustrates 
the  symbols  used  for  the  f-k  spectral  maps  shown  in  Figures  15, 

16,  17,  18,  and  19. 

Comparing  Figure  15  and  Figure  17  we  see  that  an  18% 
increase  in  the  length  of  the  vertical  array  makes  only  slight 
improvement  in  the  array  response.  An  array  with  maximum  up¬ 
hole  time  of  .55  sec  with  6  instruments  provides  an  adequate 
response  in  the  signal  band  with  the  power  to  discriminate  body 
waves  at  frequencies  above  .7  cps. 

Figures  16,  18,  and  19  show  the  array  response  for 
equally  spaced  instruments  placed  in  the  lower  2/3  of  the  well. 

The  apparent  resolving  power  is  inadequate  between  1.0  and  1.3 
cps  for  all  of  the  processors  analyzed.  The  fan  filter  appears 
to  have  the  desired  pass-band  response  with  better  dead-band 
rejection  than  the  beamed  sum  and  multichannel  deghoster. 

Several  of  the  simplest  vertical  array  processors  were 
evaluated  by  means  of  frequency-wavenumber  spectral  analysis. 

The  processors  analyzed  were  beamed  sum,  multichannel  deghost, 
and  the  fan  filter.  Adequate  responses  were  obtained  in  the 
specified  signal  pass-band  for  all  of  the  processors.  Con¬ 
siderable  differences  between  the  processors  were  observed  in 
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Figure  15.  Beamed  Sum 
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Figure  16.  Beamed  Sum 
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Figure  19.  Fan  Filter 
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TABLE  4 


Symbols  For  the  F-K  Spectral  Mapping 


DB 


Symbol 


0-1  0 

1.001  -  3  0 

6-9  6 

12-15  2 

18-21  8 

24-27 


the  dead-band,  especially  at  low  wavenumbers  corresponding 
to  surface  wave  components  in  the  noise.  The  best  dead-band 
response  was  obtained  for  the  fan  filter  using  a  spacing 
with  uphole  time  between  instruments  of  .1  seconds  from  the 
top  to  the  bottom  of  the  well. 
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H.  Hyperfine  Beamsteering  Using  a  Signal  Crosscorrela tion 
Technique 

By  beamsteering  we  mean  shifting  different  array 
channels  in  time  and  forming  a  sum.  The  beam  output  is  fre¬ 
quently  also  referred  to  as  the  "phased  sum"  or  "delay-and-sum 
output"  (Chiburis  and  Hartenberger,  1966) .  The  interchannel 
time  shifts  consist  of  those  calculated  from  an  assumed  signal 
propagation  velocity  across  the  array  and  an  assumed  direction 
of  signal  arrival  so  well  as  shifts  which  are  intended  to  com¬ 
pensate  for  variation  in  transmission  path  phase  characteris¬ 
tics  between  elements  of  the  array,  i.e.,  travel  time  anomalies. 
For  LASA,  these  time  adjustments  have  been  observed  to  depend 
on  the  signal  azimuth  and  phase  velocity,  that  is,  on  the 
signal  source  region  (Chiburis,  1966). 

In  beamsteering  the  LASA  subarrays,  only  signal  velocity 
is  used;  the  intra-subarray  travel  time  anomalies  are  negligible. 
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However,  when  beaming  subarray  beams  the  travel  time  anomalies 
are  significant:  reduction  in  signal  loss  of  approximately'  5  do 
has  been  observed  wher.  the  travel  time  anomalies  are  applied 
(Chiburis,  1966) . 

However,  signal  loss  of  approximately  1  db  has  been 
observed  in  beaming  LASA  subarrays  (Chiburis  and  Hartenberger, 

1966) ,  and  another  2-3  db  of  signal  loss  occurs  when  beaming 
the  subarray  beams.  There  appear  to  be  two  major  contributing 
factors  to  the  latter  signal  loss  which  are  variation  in  signal 
waveform  across  LASA  and  small  errors  in  the  travel  time  anomalies 

applied . 

In  this  report  we  describe  a  method  which  should  reduce 
signal  losses  due  to  the  second  of  these  causes.  This  method  can 
be  fully  automated,  provided  that  the  signal  arrival  time  is 
known  and  the  signal/noise  ratio  is  not  too  small  further  work 
wiH  be  required  to  establish  the  threshold. 

The  first  step  is  to  measure  the  differences  in  signal 
arrival  time  which  remain  in  the  subarray  beams  after  they  have 
been  time-shifted  by  both  signal  velocity  and  azimuth  and  the 
travel-time  anomalies.  Given  the  signal  arrival  time  and 
knowledge  of  the  signal  duration  based  on  prior  experience, 
two  methods  for  computing  the  signal  arrival  time  differences 
suggest  themselves. 

In  the  first  method  we  choose  one  of  the  N  subarray 
beams  as  a  reference  channel.  Compute  the  N-l  crosscorrelations 
of  the  other  subarray  beams  with  the  reference  channel,  using  a 
properly  chosen  time  window  around  the  signal  arrival  time  (at 
this  point  the  necessity  for  already  having  the  signals  almost 
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lined  up  is  apparent) .  Measure  the  lag  at  which  each  of  these 
N-l  crosscorrelations  peak.  Assume  that  these  lags  give  the 
inter-channel  time  shifts  remaining  to  be  applied. 

In  the  second  method  we  compute  the  entire  correla¬ 
tion  matrix  of  the  N  subarray  beams,  using  a  properly  chosen 
time  window  around  the  signal  arrival  time.  Measure  the  lag 
for  each  of  the  (n-l)  crosscorrelations  peak.  Assume  that 
these  represent  differences  between  inter-channel  time  shifts 
remaining  to  be  applied.  There  are  thus  ^N(N-l)  observations 
from  which  we  can  determine  the  N  time  shifts  by  least  squares, 
having  ^N(N-3)  degrees  of  freedom  (which  is  a  lot) .  Estimate 
the  standard  error  of  each  time  shift,  and  do  not  apply  any 
time  shift  which  is  not  significantly  different  from  zero. 

This  should  avoid  the  danger  of  nonsense  results  when  the 
signal/noise  ratio  is  small;  see  the  discussion  below. 

The  second  method  is  obviously  better  than  the  first 
since  it  uses  all  possible  correlations  instead  of  just  N  of 
them.  Also,  the  first  method  can  be  expected  to  fail  if  the 
reference  channel  happens  to  have  a  small  signal  or  a  signal 
waveform  profoundly  different  from  the  others  (significant 
variations  in  waveform  have  been  observed;  see,  e.g.,  Flinn 
et  al . ,  1966) .  1 

The  only  apparent  relative  drawback  to  the  second 
method  is  computing  time;  but  the  difference  in  computing 
time  between  the  first  and  the  second  method  is  actually 
quite  small.  Consider  21  subarray  beams  and  a  time  window 
150  points  long.  Then  the  first  method  requires  the  computation 
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of  twenty  lags  of  each  of  twenty  150-point  correlations,  which 
requires  a  small  fraction  of  a  second  on  the  CDC  1604-B.  The 
second  method  requires  the  computation  of  231  correlation 
functions  (our  library  routine  calculates  the  auto-correlations 
as  well  as  the  cross-correlations),  and  the  solution  of  210 
condition  equations  in  20  variables.  Forming  the  correlation 
matrix  should  require  approximately  three  seconds.  It  turns 
out  that  the  solution  of  the  least-squares  normal  equations 
for  the  time  shifts  can  be  written  down  explicitly  without  any 
matrix  multiplication  or  inversion.  Thus  the  computing  time 
is  approximately  three  seconds,  which  is  not  e;xcessive. 

Notice  that  the  problem  is  not  quite  complete  as  we 
have  stated  it  so  far:  a  constraint  is  necessary.  We  choose 
one  channel  as  a  reference  channel  and  constrain  its  time  shift 
to  be  zero;  i.e.,  we  compute  all  the  shifts  relative  to  a  refer¬ 
ence  channel  (we  emphasize  again  that  the  accuracy  of  the  least- 
squares  computed  shifts  in  no  way  depends  on  the  reference  chan¬ 
nel's  having  a  good  signal,  as  does  the  first  method). 

The  results  of  this  study  are  presented  in  the  second 
part  of  Table  5  and  in  Figures  20  through  25.  Figure  20  shows 
the  signals  from  the  three  events  used  for  computing  the  cross¬ 
correlation  matrix.  All  5  seconds  shown  were  used  to  compute 
cross-correlation  functions  sampled  as  often  as  the  data  (20 
samples  per  second)  and  out  to  a  maximum  time  lag  of  one 
second.  Figures  21,  22,  and  23  show  the  input  and  output 
traces  for  each  case;  that  is,  input  to  and  output  from  the 
vernier  beamsteering  program.  In  these  figures,  it  is  oossible 
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DEC.  9/  1967 


NOV.  H.  1965  NOV.  20,  1965 

Figure  20.  Signals  used  for  computing  the 
cross-correlation  functions 
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BEFORE  VERNIER  BEAMS fEERING  AFTER  VERNIER  BEAMSTEERING 


Figure  22.  Nov.  20,  1965  event  subarray  sums 


NOV.  H.  1965  EVENT  PPA 


Figure  24.  LASA  phased  sums  before  vernier  beams terring 


PPA  Signal  (mp) 

6.59 

13.15 

48.14 

PPA  Noise  rms  (mp) 

0.33 

0.31 

0.34 

PPA  S/N 

19.69 

42.61 

141.94 

PPAV  Signal  (mp) 

6.59 

13.57 

48.66 

PPAV  Noise  rms  (mp) 

0.33 

0.30 

0.34 

PPAV  S/N 

19.69 

45.76 

143.74 

Gain  (db) 

0.0 

0.6 

0.1 

Standard  deviation  (sec’ 

.072 

.049 

.022 

to  see  slight  time  shifts  that  were  detected  and  subsequently 
eliminated  by  the  shifting  process.  Figures  24  and  25  are 
the  sums  before  and  after  processing  respectively.  The  second 
half  of  Table  5  lists  the  numerical  results  of  the  study. 

The  last  entry  in  Table  5  is  the  standard  deviation  o* 
the  time  shifts  computed  as  described  above.  All  of  the  esti¬ 
mated  time  shifts  for  the  event  were  within  2.6  standard  de¬ 
viations  of  zero  and  were  therefore  eliminated.  The  PPAV  for 
this  event  is  thus  equal  to  the  PPA.  However,  for  two  oth  r 
events,  six  and  four  time  shifts  respectively  were  larger  than 
2.6  standard  deviations.  Application  of  these  shifts  produces 
the  differences  in  the  input  and  output  for  these  events.  At 
no  time  did  we  have  to  eliminate  a  time  shift  that  was  greater 
than  half  a  second. 

The  standard  deviation  itself  is  a  good  measure  of 
how  well  the  method  works.  The  first  event  had  the  lowest 
signal-to-noise  ratio  on  the  PPA  trace  and  the  largest  stand¬ 
ard  deviation.  The  large  standard  deviation  occurs  both  be¬ 
cause  of  the  dissimilarity  between  the  signals  from  this  event 
and  the  larger  proportionate  noise  level.  On  the  other  hand, 
another  event  had  the  highest  signal-to-noise  ratio  and  the 
lowest  standard  deviation  but  did  not  exhibit  the  largest 
gain  from  the  vernier  beamsteering  procedure.  It  appears 
that  the  anomalies  for  this  event  were  more  accurate  than 
those  for  the  last  event.  This  is  supported  by  the  lower  number 
of  non-zero  shifts. 
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These  results,  though  small  (0.6  and  0.1  db  respectively) 
are  nevertheless  considered  significant.  They  are  corrections 
to  a  refined  and  very  accurate  process,  the  phase  summation  of 
LASA  data  using  travel  time  anomalies.  Therefore  we  should 
expect  them  to  be  small.  However  what  makes  them  important  is 
that  they  result  from  eliminating  one  more  cause  of  signal  de- 
gredation  in  phased  summation.  If  one  half  a  db  can  be  saved 
by  vernier  beamsteering  only  11  subarray  outputs,  perhaps  a  full 
db  can  be  saved  by  processing  all  21.  We  have  demonstrated  here 
a  method  which  resulted  (in  one  case)  in  reducing  by  one-half 
the  signal  loss  actually  observed  (Chiburis  and  Hartenberger, 

1966) . 
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I.  Earthquake  and  Explosion  Analysis 

Standard  statistical  summary  reports  were  issued  on  the 
Nevada  Test  Site  nuclear  events,  COMMODORE,  SCOTCH,  and  KNICKER¬ 
BOCKER.  Table  6  summarizes  part  of  the  seismic  data  compiled 
from  these  analyses. 
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Other  projects  completed  included  photo  albums  for 
two  events  and  travel-time  computations  for  NTS  and  other 
sites. 

No.  of  Stations  Recording 


Event  Name 

Event  Date 

Medium 

Magnitude 

SP  Siqnals 

LP  Siqnals 

COMMODORE 

20  May  1967 

Tuff 

5.68  +  0.56 

21 

21 

SCOTCH 

23  May  1967 

Rhyolite 

5.51  +  0.75 

21 

21 

KNICKERBOCKER 

26  May  1967 

Rhyolite 

5.54  +  0.42 

16 

16 

TABLE  6.  Seismic  Data  from  Three  Nuclear  Explosions 

III.  SUPPORT  AND  SERVICE  TASKS 

A.  VELA-Uniform  Data  Services 

As  part  of  the  contract  work  statement,  the  SDL  provided 
one  or  more  of  the  following  support  and  service  functions  for 
VSC  and  other  VELA  participants: 

-  copies  of  16  and  35  mm  film. 

-  playouts  of  earthquakes  and  special  events. 

-  copies  of  composite  analog  tapes. 

-  composite  analog  tapes  of  special  events. 

-  use  of  1604  computer  for  checking  out  new  programs  or 
running  production  programs. 

-  copies  of  digital  programs. 

-  digitized  data  in  standard  formats  or  special  formats 
for  use  on  computers  other  than  the  1604. 

-  running  SDL  production  programs,  such  as  power  spectral 
density  and  array  processing  on  specified  data. 

-  digital  x-y  plots  of  power  spectra  or  digitized  data. 

-  signal  reproduction  booklets. 

-  space  for  visiting  scientists  utilizing  SDL  facilities 

to  study  data  and  exchange  information  with  SDL  personnel. 
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During  this  report  period,  54  such  projects  were 
completed  and  the  16  organizations  receiving  these  services 
are  listed  in  Appendix  A. 

B.  Data  Library 

The  data  library  contains  approximately  7400  digitized 
seismograms,  225  digital  computer  programs  and  293  composite 
analog  magnetic  tapes,  all  available  for  use  by  the  VELA-Uniform 
program. 

The  following  additions  were  made  during  this  period: 

1.  Digital  Seismograms  -  400  including 

-  data  from  two  explosions 

-  one  earthquake  recorded  at  various  stations 

2.  LASA  Data  -  9  tapes 

-  there  are  a  total  of  1396  digital  tapes  in  the  library 
including  992  field  tapes.  There  is  also  a  master 
calibration  tape  which  contains  the  magnification 
(digital  counts  per  millimicron)  of  each  sensor  for 
every  subarray.  These  magnifications  have  been  computed 
for  all  calibration  tapes  currently  in  house.  As  each 
new  calibration  is  received,  it  is  routinely  run  through 
the  new  program  CALIBR  and  added  to  the  master  tape. 

3.  Digital  Programs  -  18  including 

TFOASTRO  -  Information  processing  -  to  process  an  ASTRODATA 
Seismic  Data  acquisition  system  tape  and  output 
a  digital  to  analog  (DAC)  plot  tape. 

TFXARII  -  Data  converter  for  TFO  acquisition  system  tapes. 

This  program  makes  a  standard  library  seismogram 
tape  from  digital  tapes  in  the  format  of  TFO 
acquisition  system  tapes. 
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DEPCATE 


VFSPTRM2 


LLTOXY  - 


LOCMEV  - 


NETMEV  - 


NORFORM 


TSWITCH 


■  Intermediate  program  of  the  DEPTHMAG  LOCMEV 
sequence  designed  to  read  event  data  from  a 
DEPTHMAG  output  tape  and  prepare  an  input 
tape  for  LOCMEV. 

-  This  program  offers  an  option  to  project  a 
planar  array  of  seismometer  coordinates  onto 
the  back-azimuth  line  before  computing  the 
frequency  wave-number  power  spectra.  The 
option  of  printing  out  the  frequency  wave- 
number  power  spectra  along  a  velocity  line  is 
also  available. 

Given  the  latitude,  longitude  of  one  site,  the 
program  computes  the  x  and  y  components  in 
kilometers  from  this  site  to  other  sites  of 
known  latitude  and  longitude  within  a  seismic 
array. 

This  is  a  modification  of  LOCATECQ  designed  to 
process  events  from  an  input  tape  produced  by 
either  the  NETMEV  or  DEPCATE  programs.  LOCMEV 
computes  the  latitude-longitude  95%  confidence 
limits  for  the  input  event. 

A  modification  of  NETW0RK2  designed  to  produce 
an  output  tape  containing  event  locations  and 
selected  stations  to  be  processed  by  the  program 
LOCMEV.  The  selection  of  stations  for  input  to 
LOCMEV  are  determined  by  variable  criteria  des¬ 
ignated  on  the  input  data  deck.  LOCMEV  computes 
95%  confidence  ellipses  for  detection  of  the 
event  by  the  selected  stations. 

-  Processes  an  Astrodata  tape  containing  short 
period  data,  and  outputs  an  unpacked  library 
tape. 

-  To  prepare  a  binary  tape  for  program  LOCATE, 
containing  station  directory  and  travel  time 
and  other  tables. 
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TFOSAN  -  Processes  TF037  data,  20  samples  per  second,  less 
than  4000  points  per  channel  from  a  special  SDL 
library  tape.  This  special  SDL  library  tape  is 
made  from  processing  TF037  data  through  program 
TFSARYII,  and  then  through  program  MERGSEIS.  All 
data  to  be  evaluated  by  TFOSAN  must  follow  this 
procedure.  The  calibrations  to  be  used  as  data 
input  must  first  be  processed  by  T FOCAL. 

ARRAYPR  -  Converts  an  A-D  format  tape  into  a  UES  format 

library  tape  and  processes  the  seismograms  from 
an  array  of  seismometers. 

MAKTAPE  -  Program  searches  digitized  data  library  tapes 

(packed  and  unpacked)  for  requested  seismograms 
and  either  copies  them  (on  a  scratch  tape  in 
1604  Fortran  binary  format)  or  writes  them  in 
a  compatible  IBM  7090  format  for  use  by  other 
computer  installations. 

BMPATRN  -  To  calculate,  tabulate,  and  plot  the  azimuthal 
beam  patterns  for  an  array  steered  to  enhance 
a  particular  velocity  or  wave  number. 

LIST80  -  Fortran  63  program  to  list  columns  1-80  of 
program  or  data  cards. 

TFOCAL  -  This  program  computes  magnification  levels  in 
counts/microns  from  the  TFO  calibration  data. 

HRVKSPEC  -  This  program  computes  the  high-resolution 

wavenumber  spectrum  of  seismic  linear-array 
time  series  data.  The  spectrum  may  be  cal¬ 
culated  with  respect  to  a  particular  sensor 
in  the  array  or  it  may  be  averaged  over  all 
of  the  sensors  in  the  array. 

COMPOS IT  -  This  program  copies  the  header  record  and 

section  of  data  from  the  TFO  acquisition  system 
tape  onto  a  composite  master  tape. 
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TFSARY3  This  program  makes  a  standard  library  seismogram 
tape  from  digital  tapes  in  the  format  of  TFO 
acquisition  system,  composite  tapes  from  the 
program  COMPOSIT. 

4.  Analog  Composite  Tapes  -  4  including 

a.  Made  by  SDL 

CHARCOAL 

AKTANUM 

b.  Made  by  Geotech 

GASBUGGY 

LANPHER 


Data  Compression 


This  is  a  continuing  routine  operation,  and  production 
is  maintained  at  the  level  needed  to  meet  the  requirements  of 
the  field  operation  (LRSM  and  U.  S.  Observatories)  and  the  Seismic 
Data  Laboratory.  For  this  period  2218  tapes  were  compressed. 

D.  Automated  Bulletin  Process 

November  and  December  1967  LRSM  and  Observatory  Bulletins 
were  processed  during  this  report  period  and  forwarded  to  Geotech, 
a  Teledyne  Company,  for  checking  and  publication. 
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January  -  March  1968 

Organizations  Receiving  SDL  Data  Services 


Vitro 

Geotech 

Pennsylvania  State  University 
Colorado  School  of  Mines 
International  Seismological  Centre 
IBM 

Texas  Instruments 

General  Atronics  Corporation 

Lincoln  Laboratory,  MIT 

Southwest  Center  for  Advanced  Studies 

U.  S.  Coast  &  Geodetic  Survey 

Dallas  Seismological  Observatory 

Unitech 

Laurence  G.  Hanscom  Field 
University  oif  California 

Observatorio  San  Calixto,  La  Paz,  Bolivia 
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